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ABSTRACT 


The specific impetus for this work was a conceptual design of a 
submarine using the toroid as the pressure hull. This work is a 
continuation of the work started by Bowen (1987). As such, it is hoped 
that a better understanding of the behavior of a toroid under 
hydrostatic pressure can be realized. 

This work began with a review of efforts to solve complete toroidal 
structures. A specific toroid was then modeled in the BOSOR4 computer 
program to obtain displacements of the meridian under hydrostatic 
pressure. Functions were then derived that described the general form 
of these displacements. Using these functions as the assumed solution 
for the energy method an energy balance was made and a program was 
written to solve for the displacements of a generic toroid under 


hydrostatic loading. 
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CHAPTER 1 
INTRODUCTION 


The purpose of this work was to examine the response of a thin 
circular toroidal shell under external hydrostatic pressure. The 
approach used was to model a specific toroid in the BOSOR4 computer 
program (Bushnell, 1977) and then use the output results’ for 
displacements to generate functions that modeled these displacements. 
Using the functions obtained that model the displacements of the toroid 
under hydrostatic pressure, an analytical analysis was performed with 
these functions used as the assumed solution for the energy method. The 
energy method generated the equations necessary to solve for the 
displacements of a generic circular toroid. A program was then written 
to solve the equations generated in the energy method so_ that 
displacements for any specific toroid under hydrostatic pressure could 
be obtained. 

The analysis conducted only considers displacements in the linear 
elastic range. For the shell being analyzed the following assumptions 


were made: 


Uy The material of the shell is isotropic and homogeneous. 

2s The thickness of the shell is constant. 

S The thickness of the shell is small compared to the radii of 
curvature. 

4, The displacements are symmetric about the X-Y plane (see 
Riguneh. ll). 

oe Thin shell theory holds: normals to the undeformed surface 





remain normal. 

For a symmetrically loaded shell, with the above assumption of 
symmetric deformation, a small displacement of any point can be resolved 
into two components: v in the direction of the tangent to the meridian 
and u in the direction of the normal to the middle surface (Timoshenko 
and Woinowsky-Krieger, 1959). Therefore, the only displacements 
discussed throughout the text will be u and v. 

The toroid’s geometry is straight forward but in order to be 
consistent throughout the text the following definitions will be made 


and are illustrated in Figure 1.1: 


DEFINITIONS: 

R : Major radius of toroid. Radius of rotation about the Z axis 
of the circle of radius r to form the toroid. 

r : Minor radius of toroid. Radius of the circle which is rotated 
about the Z axis. 

h : Thickness of the shell (assumed to be constant). 

Or: Angle of rotation measured counter clockwise from the X-Y 
plane at a distance R + r from the origin. 

@? : Angle of rotation about the Z axis measured counter clockwise 


from the positive X axis. 





Figure 1.1 


Description of Geometry 








CHAPTER 2 
COMPUTER MODELING 


A computer solution for the displacements of the toroid under 
external hydrostatic pressure was sought to determine the shape of these 
displacements for a specific toroid. The computer program used to model 
the toroid was BOSOR4 (Bushnell, 1977). The BOSOR4 computer program is 
a hybrid finite difference-finite element program used to analyze 
complex shells of revolution. This program was developed by David 
Bushnell at Lockheed Missiles & Space Co., Palo Alto, California. 

Although specific reference was not found that indicated that 
BOSOR4 could adequately model a toroid, it was believed that since 
BOSORS5, a similar program that considers elastic-plastic material 
behavior, had been used to model a toroid that BOSOR4 would also 
adequately model this structure (Bushnell, 1985). 

The specific toroid that was analyzed with BOSOR4 was the toroid 
specified as the default values by Bowen (1987). The toroid geometry 
used was as follows: 

R = 8 inches 
r = 2 inches 
h = 0.02 inches 
Young’s Modulus, E = 30 x 10° psi 
Poisson’s Ratio, v= 0.3 

Modeling the toroid in BOSOR4 was accomplished by taking advantage 

of the symmetry of loading of the structure about the X-Y plane. 


Therefore, only half of the structure was discretized in the program. 





Because only half of the structure was modeled in the program, boundary 
conditions had to be established at the symmetry points. The boundary 
conditions imposed on the structure are as follows: 

@ 6 =| 0 radians 

Tangential Displacement, (v) = 0.0 
Rotation in the 6 direction, S3 = 0.0 
@ @ = TI radians 
Tangential Displacement, (v) = 0.0 
Rotation in the 6 direction, 3 =- 0.0 

The normal displacement (u) is defined as positive in the direction 
of the outward pointing normal. The tangential displacement is defined 
as positive in the direction of increasing angle @ (See Figure 2.1). 

Having the toroid modeled in BOSOR4, the program was run to 
determine the pressure required to buckle the structure. The thought 
being that if the buckling pressure predicted by BOSOR4 was consistent 
with the buckling pressures predicted by theory then the pre-buckling 
displacements generated by BOSOR4 could be used to predict the general 
form of the displacements of a generic toroid. 

The buckling pressure obtained from BOSOR4 for the toroid modeled 
differed by less than 10% from the predicted buckling pressure as 
presented by Sobel and Fltigge (1967). Since the buckling pressure 
obtained was consistent with the predicted pressure and the predicted 
buckling pressure was substantiated experimentally by Almroth, Sobel and 


Hunter (1969), it was concluded that the results from BOSOR4 for the 


10 





Eigures 2. . 


BOSOR4 Modeling 
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pre-buckling displacements could be used to predict the general form of 
the displacements of a generic toroid loaded under uniform hydrostatic 
pressure. 

In any finite difference or finite element program an investigation 
into the number and spacing of the nodes used to describe the structure 
must be undertaken to ensure that the structure is adequately described. 

For this investigation, equal spacing of the nodes was used throughout. 

To investigate the effect of node distribution, the number of nodes 
were doubled and thus the spacing was reduced by a factor of two. The 
results for both the magnitude of the buckling pressure and the shape 
and magnitude of the displacements did not change with the increase in 
the number of nodes. It was therefore concluded that the structure was 
adequately described and that the results obtained were valid. 

Figure 2.2 is a graphical presentation of the normal displacement 
(u) and figure 2.3 is a graphical presentation of the tangential 


displacement (v) obtained from the BOSOR4 output for the toroid modeled. 


Wy 





u CGdnohes) 


Figure 2.2 
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CHAPTER 3 
CURVE FITTING 

The next step in this study was to generate a function that would 
analytically describe the normal and tangential displacements of the 
toroid under hydrostatic loading. 

First a check was made, using the BOSOR4 computer program, to 
ensure the displacements were linearly dependent on the applied 
pressure, as should be the case since this study was only considering 
behavior of the toroid in the linear elastic region. As can be seen 
from figure 3.1 the normal displacement is indeed linearly dependent on 
the applied pressure. Although not included here, the tangential 
displacement is also linearly dependent on the applied pressure. 

To generate a function that would adequately describe the 
displacements, a form of the displacements must be assumed and then a 
curve fitting technique must be used to fit the data to this assumed 
form. 

since the displacements are considered to be symmetric about the 
X-Y plane and considering the way in which the positive direction for 
the normal displacement (u) and the tangential displacement (v) were 
defined in Chapter 2, the following statements can be made: 

For the normal displacement (u), 

u(@) = u(-6) 3a 1) 
which is the definition of an odd function. 

For the tangential displacement (v), 


v(@) = -v(-6) (322) 


Ge, 





u Gnohea) 


Figure 3.1 
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which is the definition of an even function. 

The functions assumed to describe the displacements were Fourier 
Series for both the normal (u) and tangential (v) displacements. Since 
both the normal and tangential displacements have a period of 2xr (the 
arc length of the circular cross section) and noting the above 
statements about even and odd functions, the Fourier Series for the 


displacements can be written as follows: 


fs @] 


u(é) = = + % an cos [= eae oe 
n=1 
v(@) = JY bn sin [Pe (3.4). 
n=1 


For the displacements to be symmetric about the X-Y plane, or a 
half period of the function, x and L can be defined as follows: 
x= ré , 


L = xr 


From these definitions, u and v can be rewritten as: 


u(@) = = + y an cos (né) Cares) 
n=1 
v(6) = }) bn sin (né) (OUD 
n=1 


While the Fourier series described above will give an exact 
representation of the displacements, the goal of a curve fitting routine 


should be to determine how many terms of the series must be used to 


7 





adequately describe the function. To determine the number of terms 
required to describe the displacements a least squares curve fitting 
program was written to be used in conjunction with the output data from 
BOSORG. The programs were adapted from a curve fitting routine 
presented by James, Smith and Wolford (1977) and are included in 
Appendix A. 

Using these programs to determine the coefficients of the Fourier 
Series it was found that five terms were inadequate to describe the 
displacements, as can be seen graphically in Figure 3.2, and it was 
determined that ten terms were required to adequately describe the 
normal displacements, (See Figure 3.3), and nine terms were required to 
describe the tangential displacements. Although only the normal 
displacement (u) is presented in Figures 3.2 and 3.3 similar results 
were obtained for the tangential displacement (v). 

With the number of terms determined that are required in the 


Fourier series u and v can be written in final form as: 


9 


u(@) = = + ), an cos (né) OS ya 
n=l 
9 
v(6) = } bn sin (né) Coe 
n=l 


These definitions of u and v will be used in the next section as the 


assumed forms of the displacements of the toroid. 
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CHAPTER 4 
ENERGY METHOD 
4.1 Introduction: 

To determine the displacements of a structure when it is subjected 
to an external load, the external load must be balanced by the internal 
forces of the structure caused by the external load. 

For this study, the method of minimizing the total potential 
energy, as described by Shames and Dym (1985), was used. The energy 
method was chosen in hopes of avoiding the singularities in the solution 
at @=t ,(Reissner, 1963 and Senjanovic , 1972). The total potential 
energy is defined as: 

T=-U- W Cave) 
where U and W are defined as; 
U : internal energy of the structure, 
W : potential energy of the applied loads. 
To determine the minimum potential energy, and thereby determine the 
displacements caused by the external load, the first variation of the 


total potential energy is set to zero. 


(1) 


él] = 0 (41,2) 


The first variation of the total potential energy can be defined as the 


following partial differential equation: 
= =0 (4.1.3) 


where x is any general parameter used to describe U and W. 
1 


Since the first variation of the total potential energy is set to zero, 


2a: 





the equation to be solved to determine the displacements of the 


structure can be defined as; 


Ul ge Uh Ae 


In order to solve the above equation for the displacements caused by the 
applied load, U and W must be defined as functions of the displacements. 
As described by Bushnell (1984), the strain energy of the structure 


can be defined as follows: 


eS eee C4 Ao) 
where; 
ue : Strain energy due to elongation and changes in 
curvature, 
U_: strain energy due to constraints or boundary conditions. 


Cc 


For the toroid, the symmetry about the X-Y Plane will be taken 
advantage of and only half of the structure will be analyzed. The 
following additional geometric definitions are required for the analysis 


(See Figure 4.1). 


DEFINITIONS: 
a: Ratio of the major and the minor radius of the toroid. a= = 
r_: Distance from the Z axis to the meridian of the toroid. 


r. - R+ rcos(é) = ell a+ cos (6)) 


r : Radius of curvature in the @ direction. 


ae R + r cos(6) 


~. ~ Gos(#)  cos(d) 


2 





Figure 4.1 
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4.2 Strain Energy: 
The strain energy due to elongation and changes in curvature can be 
defined as: 


i ; 
oo aaa ec wmey = isd 5 o Cory 


The strain energy due to constraints, ue , will be addressed later. 

Because this study is only considering hydrostatic loading of the 
toroid, the loading will always be normal to the surface of the toroid. 
Therefore, there will not be any shear stresses or shear resultants from 
the loading imposed. 


Caecum = 0 4 tt 4 eG), 
ij ij 


The @ and ¢ directions are thus principle directions and the strains can 
be represented by reference surface strains and changes in curvature. 
The strain energy can be defined in terms of forces and moments which 
are defined as an integral of the stresses through the thickness of the 
shell. The strain energy can then be represented by a surface integral 


as follows: 


1 
Ue I, (N,, per iets. 1F\. ) ds Gee?) .. 


For the case of hydrostatic loading the strain energy becomes: 
1 
U=5 f(g, ~ Ngfg J + (Mon, + Mog ) Jas (a3) 
where the forces and moments are defined as follows (Timoshenko and 


Woinowsky - Krieger, 1959): 


Care 
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M,-|,e,2 (1-2) & (4.2.6), 


z 
My = |, pe (1-E je (4.2.7). 


Since the assumption has been made that normals to the undeformed 
middle surface remain normal, the strains as a function of thickness can 
be expressed in terms of reference surface strains and changes in 


curvature (Timoshenko and Woinowsky-Krieger,1959). 


fp “ €p - Zp (4.2.8) 
t 
an - Z@ aS) 
a ee: as 
The terms, € g and 7 , represent the total strain in the 6 and ¢ 
t t 


directions and «, and «, represent the changes in curvature of the 


y ¢ 


reference surface. 


From Hook’s Law the stresses can be defined as follows: 


pa tes 
¢, = €, + ve (4.2.10) 
: (1-v*) on a 
¢,-= ee €, + veE Cie? ais) 
¢ Que) a ye 
substituting the expressions for € and 27 ; 
t t 
og-— [etry -2( ny +e, ) | (Gr2em2):. 
Gisus) 
v4 = + 2 + VE, -z( we + UK ) Ci 23) . 
Cie) 


Taking the middle surface as the reference surface and neglecting 


Z Z : : ; 
the terms, = and > as small as compared to unity, integration over 
2 


25 





the thickness of the shell can be represented by: 

















i E ; 
N, = €, + vé, -Z| K, + VK | dz (4.2.14), 
q -h/2 45m q ¢ j ¢ 
Eh 
N. = —— €, + vE Ca Zo). 
6 Glens 6 v) 
+h/2 5 ; 
N, = €, + vé, -Z] K, + VK dz (On 10) 
? -h72 (1-v*) , q ° ” 
Eh 
N, = €e, + ve CA ee? ).., 
¢ (1-v?) ¢ Y 
+h/2 E ; 
M. = . + ve, -Z|K, + VK Zaz Cae Se 
g -h7/2 Gleue) @ ¢ Y ¢ 
3 
M = ae Kg + wo, G42 92 
ace) 
—_— 
M, = . + ve, -Z|K, + VK Zaz Cia PAS 
? meyer AC Leu) ? ‘ ? : 
3 
ae | 
qM, = os K, + UK CGn2 2): 
? Ia@leme yb : 


Combining the expressions for the forces and the moments 
substituting these back into the expression for the strain energy, 


strain energy becomes: 





U = ie Bn ( AS + 2ve.€, + oe ) - 
$s Z I. Gleus) 6 6 p v) 


Eh° 
eles a 
4.2.1 Strains: 


The next step is to define the surface strains in terms of 


displacements u and v. This will be accomplished by following 
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2 2 
( Kp + 2UK,.K, + S ) ds CLE a 


and 


the 


the 


the 





presentation by Timoshenko and Woinowsky-Krieger, (1959). 
Considering an element AB of the meridian, (See Figure 4.2), the 
strain in the @ direction can be defined as follows: 


é, due to v, 


8 
Ov 
Vim Y 1 av Us, Ga 
: r dé r 86 eam 


Ey due to u can be represented by the average change in length of 


the element divided by the original length, 


1 du 


. Stra ve vee 
ary ea ellie) 
The component, : = dé , will be neglected as a small quantity of higher 


order. The total strain in the 6 direction can then be represented as: 


Ih Ov 
oe u + Ey] | (Ame 13). 


The strain in the ¢ direction can be defined as the change in i 
due to u and v divided by the original length of the element (See Figure 
4.3). 


A r. = ( u cos(@) - v sin(@) ) d@ , and 


( ucos(@) - v sin(6) } d¢ 
nS ee OAL A 
é ro 
1 


“¢ ~ ¥ (a + cos(6)) 


fu cos(8) - v sin(#) | (4.22 eles) 2 
4.2.2 Change in Curvature: 
The change in curvature of the middle surface will now be defined. 
Again considering the element AB of the meridian, (See Figure 4.2), 


the rotation in the @ direction can be defined as follows: 
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Figure 4.2 
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Figure 4.3 


Radial Displacement 
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p-— 4 Cp lnral) 
peaue EO. Ww, 
du 
- = ge 
p-—2 = Bag CE ae 


Mae total rotation in the @ direction is; 


1 du 
B, _ = [ : se | ag (Agee? 5.) 


The change in curvature in the @ direction can be represented by the 
change in rotation in the @ direction divided by the undeformed arc 


length in the @ direction. 


6 2 
r 


op 2 
8 iL Ov Qu 
a ae] (4.2.2.4) 


x 
i 
rire 


Because the assumption has been made that the deformations will be 
symmetric, the rotation in the @¢ direction will be opposite to the 
rotation in the @ direction in magnitude and rotated into the ¢ surface 


by sin (@). 


By = - By sin (6) = | = - Vv sin (6) d¢ (an 2 25) 


Although the rotation in the ¢ direction is constant in ¢@ , it does vary 
in the @ direction. Therefore, the change in curvature due to loading 
can be represented by the rotation in the ¢ direction divided by the 


undeformed arc length in the ¢ direction. 


Pg 
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and substituting the expression for a 


Ce aaa 


sin(6) du 
=e (8+) 


r*(a + cos(6)) 


4.2.3 Constraint Energy: 

A complete toroid has no constraint or boundary conditions when 
loaded only by hydrostatic pressure. However, since only half of the 
toroid is being analyzed here, boundary conditions must be imposed at 
the symmetry points, 6 = 0 and 6 = x. The boundary conditions that are 
required at these two points are: 


The vertical displacement is zero, v = 0, 


and the rotation in the @ direction is zero, = =- Q. 


As Bushnell (1984) points out the constraint energy can be accounted for 
with the introduction of Lagrange multipliers and can be represented as 


TOLlows: 
U =i eo + X v(O0) + A oe +X vin) C427 Se) 
c Oud 8 Ov ud 8 Av ne es : 


In this formulation of the constraint energy the Lagrange multipliers 
are unknowns and are solved for along with the displacements when the 
total potential energy is minimized. The actual form of the constraint 


energy will be discussed later. 


4.3 Potential Energy of Applied Loads: 
The potential energy of the applied loads is simply the work done 


by these loads, which can be represented by the applied load multiplied 


Sik 





by the deflection in the direction of the applied load. For the case of 
hydrostatic loading, the applied load is always normal to the surface 
and the displacement in the direction of this load is -u for this study. 
Therefore, the potential energy of the applied hydrostatic load can be 


represented as: 
W=- ie Eee ds Camo) 


where P is the external hydrostatic pressure applied and is 


positive in the direction of -u. 


4.4 Minimization of the Total Potential Energy: 
From section 4.1 it was shown that to minimize the total potential 


energy the following equation needed to be solved; 
a re Eis) 


and it was shown that the above equation could be represented by 


the following equation: 
ees (= (4.4.2). 


To present the form of the above equations the following simplifications 


in notation will be used: 


ee ces + «ee 
ag ° 397 ag 
3 
aa a0 pe —Eh e, 
(l-v') 2 Olaus } 
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4.4.1 Internal Energy: 


Since it has been shown that none of the terms in the internal 


energy depend on ¢ the surface integral can be rewritten and integrated 


with respect to ¢. 
x 26 
1 Eh 2 2 
U =-- | (« + 2vée.€, + € - 
8 2 | Glew § 6 g 4) 


Eh°> 2 


2 
pa.08) 8 + 2unge, + 3) eC a + cos(@)}d¢ rdé Cea) 





Which when integrated with respect to ¢ becomes, 





2 - Eh 2 2 

U =a Xr | (« + 2ve .€, + € } - 

s ee 6 6 ¢ g 
Eh° 


(x, + 2UK ok + 3] ( at cos()} dé CARA 2m 
¢ 
1 YGlneoe) 


Using the above simplifications in notation and substituting the 


expressions presented in section 4.2.1 and 4.2.2 the strain energy 


becomes: 


7 
U =" r? | 2 ( Cae 2uv'’ + vi? } + 
s r 


Tae ( (uP tuv’ ye0s(e) = (uv + wv! sin(#)) 


Sy oseee ( u’cos*(6) - 2uv cos(@)sin(6)+ 


v’sin*(6) ) - 2 (v'- 2 oe ne ay + 


2v sin(6) heel ' feat oe 
(a + cos(@)) ( a a es a 


a3 





sin*(@) 


me ale bere 7}}] (a+cos(4)) 30 (CE A Aes oe 


The above equation represents the strain energy, ee However, for 


minimization of the total potential energy the expression that is needed 





fe 
Ox 
ou n 
s du dv’ du yeh 
a n | [ ¢ ( Max, + as +v' ax? ey Ax ) + 
at z i 
2v dv' du GY yeu 
ae tay (C2 Nex, +v' 5z 008(8)- (ugg RE aire +V ax )sin(6)] 
Ih dv 
a — a — 
(Gas Saat) a cos *(6)- 2(ug Fe ae ~ )cos(#)sin()+2vz ax, 5+" *a))] 
- 2 (av - acy 3 tu’! ~ "+ Qu sar ) + 
wes i ar i 
_2 v sin(é) vise HOMMEMOR, 27 OMe OU We, OU | OU 1 BV 
(a + cos(6))\” ax" ax, "OX, V ax. os ax. ie ax ax oi ad ij 
Z 
sin’ (4) ow 8Y 4484") 9 3¥ 
C5 CIB ee ax 2(u Bx. *¥5x )*2¥%Gy)] (ete0sc09) 20 (440.1 4) 


The above equation represents the first term on the right hand side of 
the following equation that will be needed to minimize the total 
potential energy and thereby solve for the displacements, 


3W dU. dU 


¢c 
pan oe oe (4.4.1.5). 
i i i 


4.4.2 External Energy: 
Similar to the internal energy, the external energy does not depend 
on @ and the surface integral can be rewritten and integrated with 


respect to ¢. 


W=- | | Pou r( a + cos(é) } dé r dé (GXis2, 1) 
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Which when integrated with respect to ¢ becomes, 
> 
Wom -2Par | u (a+ cos(é) ) dé (Grae ee))e 


For minimization of the total potential energy the expression that is 


needed is ok 


ox 
i 


OW 2" du 
ax, -2Pxrr | ax, | a+ cos(@) ) dé Cae a 23) 
This completes the expressions needed for the minimization of the total 


potential energy, with the exception of the constraint energy, U 
Cc 


4.5 Assumed Functions: 

In Chapter 3 it was shown that the displacements u and v could be 
represented by a Fourier series. These functions of @ will be used as 
the assumed form of the displacements and the coefficients of the 
Fourier series will be solved for from the equations presented in the 
minimization of the total potential energy. The assumed functions for u 
and v are different from those presented in Chapter 3 in that the 
constant term in the expression for u has been included inside the 


summation. 


9 


u(é) = ), an cos (né) Cole oan, 
n=0 , 
9 
v(6) = } bn sin (n6) Cee) 
n=0 


The above expressions for u and v involve 19 unknown coefficients that 


will have to be solved for to describe the displacements. 


? 


Using these expressions for u and v the constraint energy, U 
Cc 
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discussed in section 4.2.3 can now be addressed. The constraints necded 
were the boundary conditions at the symmetry points, @ = 0 and 6 = nm. 
These boundary conditions are: 
v= 0. Since sin(né) = 0, for 6 = O and sin(nz) = O, for all no, 
these boundary conditions are satisfied by the assumed form of v 
and additional constraints are not needed. 


9 
du = 0. Since aks = - yn an sin (né), then gu = 0 at 6 = O and 


06 a6 dé 
n#0 
6 = x for the same reasons as stated above for v. 
Since no additional constraints are needed to meet the boundary 


conditions, Lagrange multipliers are not needed and the constraint 


energy, UL - 0. 


4.6 Solution: 

With the assumed forms of the displacements defined, the problem is 
now restricted to solving for the 19 unknown coefficients of the Fourier 
series. Using the following substitutions the equations can be set up 


to solve for these coefficients: 


9 


u= ) an cos (né6) (4.6.1), 
n=0 
9 
ul = - yan an sin (né6) Case 2 )e 
n#=0 
9 
u'’ = - yn? an cos (né) Cayce). 
n#0 
9 
v= ) bn sin (né) (4.6.4), 
n=0 
9 
vi = yn bn cos (né) (426 25). 
n#0 
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Substituting the above expressions into the following equation will then 
represent the 19 equations needed to solve for the 19 unknown 
eoefficients: 


dW gk 


= GROG a, 
i i 


Which can be represented in matrix notation as follows; 


omle= lr) [x ] (4.6.7) 


where the above terms are defined as: 


D: a419x 1 column consisting of the terms, a 
i 
aU 
F ;:; a19x 19 matrix consisting of the terms, cee 
i 
x > a1l9x 1 column consisting of the terms, an and bn. Where x, 


through xX, are the 10 an’s and x through xX, are the 9 
bn’s. 
To solve for the X'S, the coefficients of the Fourier series, the 


following matrix operation needs to be performed: 


fs ee) Camere) 


The individual elements in F can be determined by performing the 


dU 

integration from 0 to x of the expression for Aa as presented in 
i 

section 4.4.1. Similarly , the individual elements in D can be 


: : : : aw 
determined by performing the integration of the expression for 5x «as 
i 
presented in section 4.4.2. To perform the above integrations, as well 
as inverting the F matrix and solving for the x ’s, a program was 
i 


written. This "Toroid Program" is listed in Appendix B. 


oF 





The "Toroid Program" was run to compare the results with those 
obtained using BOSOR4. As can be seen from Figures 4.4 and 4.5, the 
results obtained for the displacements are comparable to those obtained 
from BOSOR4 and are of the same general shape. Additionally, when the 
internal energy of the results are compared, the "Toroid Program" is 
only 4.3% greater than that of the BOSOR4 solution. This difference can 
be accounted for by the selection of the functions chosen to represent u 
and v. The functions chosen for u and v only model the displacements 
and do not model the change in displacements or the curvature of the 


toroid. 


4.7 Forces and Moments: 


The next step was to look at the predictions for the forces, N, and 


i] 
Ms , and the moments, Mi and my , and compare them with the results from 
BOSORG. 
The forces and moments predicted by the "Toroid Program" did not 
corollate well with the results obtained from BOSOR4. As can be seen 


from Figures 4.6 and 4.7, the average of the predicted forces is 
approximately the same as the forces obtained from BOSORS. However, as 
can be seen from Figures 4.8 through 4.11, the predicted moments were 


several orders of magnitude less than those obtained from BOSOR4. 
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NORMAL DISPLACEMENT 
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Figure 4.4 
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TANGENTIAL DISPLACEMENT 


Figure 4.5 


Tangential Displacement (v) 
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Figure 4.6 


Ny vs 6 
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Figure 4.7 
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Figure 4.8 
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MERIDIAN MOMENT (lb—in) 


Figure 4.9 


M, vs @ (Toroid Program) 
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CIRCUMF. MOMENT (lb-—in) 


Figure 4.10 
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Figure 4.11 


a vs §@ (Toroid Program) 
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CHAPTER 5 
CONCLUSIONS 

As was shown in Chapter 4, using the assumed functions for the 
displacements in the energy method produced results for. the 
displacements that were consistent with those obtained from BOSOR4. 
However, it was also shown that the "Toroid Program" was very poor at 
predicting the forces and moments of the toroid. 

The differences can be explained by looking at what the assumed 
functions used for the displacements represents. The assumed functions 
for u and v are fairly good models of the displacements but were not 
modeled to meet the slope or curvature of the displacements except at 
the end points, where 6 = 0 and 6 = ff. This resulted in good 
predictions for u and v but inaccurate predictions for u’, u’’ and v’, 
all of which are contained in the expressions for the forces and 
moments. 

Some thought has been given to what other assumed functions might 
be used to model the displacements that would more accurately predict 
the forces and moments. It is believed that a complete Fourier series 
representation, instead of just an odd or even series, might better 
predict the slope and curvature of the displacements and therefore the 
forces and moments. To ensure that the complete Fourier series meets 
the boundary conditions, Lagrange multipliers would be required to 
describe the constraint energy (See section 4.2.3). This is an area 


that is left for further investigation. If a complete Fourier Series 


47 





representation does not adequately predict the forces and moments the 
next step would be to write some kind of finite element or finite 
difference program to solve this problen. 

Since the program written does model the displacements and produces 
results that are consistent with those obtained from BOSOR4, the "Toroid 
Program" can be used as a preliminary design tool without having to 
resort to a very complex and expensive computer code such as BOSOR4. 
Unfortunately the "Toroid Program" is restricted to analyzing only 
toroids under hydrostatic loading. Complexities in describing the 
geometry of a general shell of revolution will quickly lead the designer 
to the more complex computer codes like BOSOR4. 

It is also believed that the program written could be modified to 
incorporate nonlinear terms and then be used to predict buckling loads. 


This is again an area that will be left for further investigation. 
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APPENDIX A. 


CURVE FITTING PROGRAMS 


The programs listed were used in conjunction with the out put from 
BOSOR4 to determine the coefficients of the assumed Fourier Series. 
Both programs are written in Fortran 77. 

Program 1 determines ten coefficients for an odd Fourier Series 
that is symmetric about @ equal to m”, to fit the data from BOSOR4 for 
the normal displacement of the shell meridian. 

Program 2 determines nine coefficients for an even Fourier Series 
that is anti-symmetric about @ equal to ”, to fit the data from BOSOR4S 


for the tangential displacement of the shell meridian. 


PROGRAM 1 LISTING: 


C CURVE FITTING PROGRAM 
C LEAST SQUARES CURVE FITTING 
C234567 
PVEEGER I,J,M,N 
Rae erOoe lou) 7 ri 100,10) ,FT (10,100) ,Ac10,11) ,BC(10) 
Renee ho nto, F5,F6,F/,F8,F9,F10,C(10) 
Ebner ro, FS, FS ,F6,F/ ,F8,F9,F 10 
N = 100 
M = 10 
C READ X-Y VALUES OF DATA POINTS 
DO 10 I=1,N 
READ (UNIT = 2,FMT = *) X(I) 
READE UNIT ge, FMT = *) Y(1) 


30 








10 CONTINUE 
C GENERATE THE F MATRIX 
DO 20 I=1,N 
Ber) = FCC! )) 
F(I,2) = F2(X(I)) 
CN 1S GO) 
F(I,4) = F4(X(I)) 
F(I,5) = F5(X(I)) 
F(I,6) = F6(X(I)) 
F(I,7) = F7(X(1)) 
F(I,8) = F8(X(I)) 
F(I,9) = F9(X(I)) 
F(I,10) = F10(X(I)) 
20 CONTINUE 
C GENERATE THE TRANSPOSE OF THE F MATRIX 
DO 30 I=1,N 
DO 30 J=1,M 
FT(J,1) = F(I,J) 
30 CONTINUE 
C DETERMINE COEFFICIENT MATRIX A OF SIMULTANEOUS 
C EQUATION SYSTEM 
CALL MATMPY(FT,F,A,M,N,M) 
C DETERMINE COLUMN OF CONSTANTS FOR SIMULTANEOUS 
C EQUATION SYSTEM 
CALL MATMPY(FT,Y,B,M,N,1) 
DO 40 I=1,M 
A(I,M+1) = B(I) 
40 CONTINUE 
C DETERMINE A(n) VALUES BY SOLVING SIMULTANEOUS EQUATIONS 
C USING CHOLESKY METHOD 
MPl =M+1 
CALL CHLSKY(A,M,MP1,C) 
C WRITE OUT THE A(n) VALUES 


onl 





C 


WRineeCUNLE = 4, FMT = *)(1,C(1),I=1,M) 


END 


C DETERMINES MATRIX C AS PRODUCT OF A AND B MATRICES 
SUBROUTINE MATMPY(A,B,C,M,N,L) 
REAL A(M,N),B(N,L),C(M,L) 


10 


DO 10 I=1,M 
DO 10 J=1,L 
C(1,J).= 0. 
DO 10 K=1,N 


C(I,J) = C(1I,J)+ACI,K)*B(K,J) 


CONTINUE 
END 


DEFINE THE FUNCTIONS 


DEFINE THE A(0) FUNCTION 


REAL FUNCTION F1(X) 
REAL X 

Fl = 1.0 

END 


DEFINE THE A(1) FUNCTION 


REAL FUNCTION F2(X) 
REAL X 

F2 = COS(X) 

END 


DEFINE THE A(2) FUNCTION 


REAL FUNCTION F3(X) 
REAL X 

F3 = COS(2.*X) 

END 


DEFINE THE A(3) FUNCTION 


REAL FUNCTION F4(X) 
REAL X 


Sf 





F4 = COS(3.*X) 
END 

DEFINE THE A(4) FUNCTION 
REAL FUNCTION F5(X) 
REAL X 
F5 = COS(4.*X) 
END 

DEFINE THE A(5) FUNCTION 
REAL FUNCTION F6(X) 
REAL X | 
F6 = COS(5.*X) 
END 

DEFINE THE A(6) FUNCTION 
REAL FUNCTION F7(X) 
REAL X 
F7 = COS(6.*X) 
END 

DEFINE THE A(7) FUNCTION 
REAL FUNCTION F8(X) 
REAL X 
F8 = COS(7.*X) 
END 

DEFINE THE A(8) FUNCTION 
REAL FUNCTION F9(X) 
REAL X 
F9 = COS(8.*X) 
END 

DEFINE THE A(9) FUNCTION 
REAL FUNCTION F10(X) 
REAL X 
F10 = COS(9.*X) 
END 


33) 





SUBROUTINE CHLSKY(A,N,M,X) 
REAL A(10,11),X(10) 
C CALCULATE FIRST ROW OF UPPER UNIT TRIANGULAR MATRIX 
DO 10 J=2,M 
A(1,J) = A(1,J)/A(1,1) 
10 CONTINUE 
C CALCULATE OTHER ELEMENTS OF U AND L MATRICES 
DO 60 I=2,N 
J =I 
DO 30 II=J,N 
SUM = 0. 
JM1 = J-l 
DO 20 K=1,JM1 
SUM = SUM+A(CII,K)*A(K,J) 
20 CONTINUE 
Aree) = ACI1,J)-SUM 
30 CONTINUE 
IPl = I+l 
DO 50 JJ=IP1,M 
SUM = 0. 
IM1 = I-1l 
DO 40 K=1,IM1l 
SUM = SUM+A(1,K)*A(K,JJ) 
40 CONTINUE 
A(I,JJ) = (A(I,JJ)-SUM)/A(I,I) 
50 CONTINUE 
60 CONTINUE 
C SOLVE FOR X(I) BY BACK SUBSTITUTION 
X(N) = A(N,N+1) 
L = N-l 
DO 80 NN=1,L 
SUM = 0. 
I = N-NN 
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70 


80 


IPl = I+] 

DO 70 J=IP1,N 

SUM = SUM+A(1I,J)*X(J) 
CONTINUE 

X(I) = A(I,M)-SUM 
CONTINUE 

END 


2)3) 





PROGRAM 2 LISTING: 





C CURVE FITTING 
C LEAST SQUARES CURVE FITTING 
C234567 


INTECER] [J M,N 

Rpm LOO} ar (LOO) \F(1L00,9) > FT(9, 100) ,AC9, 10) ,BC9) 
Rewer, Fo, F4,F5,FO,F),F8,F9,G(9) 

BOGERNAL Fl, F2,63, Po FS, Fo,F/,F8,F9 

N = 100 

M=9 


C READ X-Y VALUES OF DATA POINTS 


10 


DO 10 I=1,N 
READ (UNIT = 2,FMT = *) X(I) 
READ (UNIT = 5;FMT = *) Y(1) 
CONTINUE 


C GENERATE THE F MATRIX 


20 


DO 20 I=#1,N 
Gre ee FL CXC.) ) 
FG =F 2 (x(t) 
Peo ye— FICK(L)) 
DCist) = FSCXC1)) 
Gro). = Fo(x(1)) 
FtL,6) = F6(X(1)) 
BC? ) = F7 (X(1)> 
Foro) = F8(X(1)) 
PCL, 9) = F9CX(1)) 

CONTINUE 


C GENERATE THE TRANSPOSE OF THE F MATRIX 


30 


DO 30 I=1,N 
DO 30 J=1,M 

Pie =F 1,J) 
CONTINUE 
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40 


10 


DETERMINE COEFFICIENT MATRIX A OF SIMULTANEOUS 
EQUATION SYSTEM 

CALL MATMPY(FT,F,A,M,N,M) 
DETERMINE COLUMN OF CONSTANTS FOR SIMULTANEOUS 
EQUATION SYSTEM 

CALL MATMPY(FT,Y,B,M,N,1) 

DO 40 I=1,M 

AG Mel) = Bl) 

CONTINUE 
DETERMINE B(n) VALUES BY SOLVING SIMULTANEOUS EQUATIONS 
USING CHOLESKY METHOD 

MPl=M+1 

CALL CHLSKY(A,M,MP1,C) 
WRITE OUT THE B(n) VALUES 

WRITE (UNIT = 4,FMT = *)(1I,C(I),I=1,M) 

END 


DETERMINES MATRIX C AS PRODUCT OF A AND B MATRICES 
SUBROUTINE MATMPY(A,B,C,M,N,L) 
REAL A(M,N),B(N,L),C(M,L) 


DO 10 I=1,M 

DO 10 J=1,L 

C(I,J) = OQ. 

DO 10 K=1,N 

GC(l-J) = C(1,J)+A(1,K)*B(K,J) 
CONTINUE 

END 


DEFINE THE FUNCTIONS 
DEFINE THE B(1) FUNCTION 
REAL FUNCTION F1(X) 


REAL X 


5/ 








Fl = SIN(X) 
END 

DEFINE THE B(2) FUNCTION 
REAL FUNCTION F2(X) 
REAL X 
F2 = SIN(2.*X) 
END 

DEFINE THE B(3) FUNCTION 
REAL FUNCTION F3(X) 
REAL X 
F3 = SIN(3.*X) 
END 

DEFINE THE B(4) FUNCTION 
REAL FUNCTION F4(X) 
REAL X 
F4 = SIN(4.*X) 
END 

DEFINE THE B(5) FUNCTION 
REAL FUNCTION F5(X) 
REAL X 
F5 = SIN(5.*X) 
END 

DEFINE THE B(6) FUNCTION 
REAL FUNCTION F6(X) 
REAL X 
F6 = SIN(6.*X) 
END 

DEFINE THE B(7) FUNCTION 
REAL FUNCTION F7(X) 
REAL X 
F7 = SIN(7.*X) 
END 

DEFINE THE B(8) FUNCTION 
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REAL FUNCTION F8(X) 
REAL X 

F8 = SIN(8.*X) 

END 


C DEFINE THE B(9) FUNCTION 


REAL FUNCTION F9(X) 
REAL X 

F9 = SIN(9.*X) 

END 


SUBROUTINE CHLSKY(A,N,M,X) 
REAL A(9,10) ,X(9) 


C CALCULATE FIRST ROW OF UPPER UNIT TRIANGULAR MATRIX 


10 


DO 10 J=2,M 
A(1,J) = A(1,J)/A(1,1) 
CONTINUE 


C CALCULATE OTHER ELEMENTS OF U AND L MATRICES 


20 


30 


DO 60 I=2,N 

J= I 

DO 30 II=#J,N 

SUM = 0. 

JM1 = J-1 

DO 20 K=1,JM1 

SUM = SUM+A(II,K)*A(K,J) 
CONTINUE 

A(II,J) = ACII,J)-SUM 
CONTINUE 

Ved 1 

DOes0 J J=1P1 M 

SUM = 0. 

IM1 = I-1l 

DO 40 K=1,IM1 

SUM = SUM+A(1,K)*A(K, JJ) 


ao 





40 


50 


60 
C 


70 


80 


CONTINUE 
A(I,JJ) = (ACI,JJ)-SUM)/A(TI,T) 
CONTINUE 
CONTINUE 


SOLVE FOR X(I) BY BACK SUBSTITUTION 


X(N) = A(N,N+1) 
L=N-l 

DO 80 NN=1,L 

SUM = 0. 

I = N-NN 

Jieihe NE ok 

DO 70 J=IP1,N 

SUM = SUM+A(I,J)¥*X(J) 
CONTINUE 

X(I) = ACI,M)-SUM 
CONTINUE 

END 
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APPENDIX B. 


TOROID PROGRAM 


This program is set up to determine the ten coefficients of the 
even Fourier Series for the normal displacement, (u), of the toroid and 
the nine coefficients of the odd Fourier Series for the tangential 
displacement, {v}, of the toroid. 

To determine the shape and magnitude of the normal and tangential 


displacement of the toroid the following equations are used: 


g 
u = ) a cos(né) 
n=0 i 


9 
v=) b_sin(né) 


n=l 


The program was written in Fortran 77. The subroutines used were 


adopted from those presented by Dyck, Lawson and Smith (1984). 


PROGRAM LISTING: 
C TOROID PROGRAM 
INTEGER I,N 
Remi on id), D( 19) ,C(19) 
REAL E,V,RB,RL,H,P,Al1,A2,A3,A 
LOGICAL ERROR 
EXTERNAL FA,FA1,FB,FB1,SIMP 
PRINT *,’This Program will calculate the coefficients’ 


PRINT *,’of the Fourier Series for the Normal and’ 
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PRINT *,’Tangential displacements of a thin, circular’ 
PRINT *,’toroidal shell when loaded by uniform’ 
PRINT *, ‘external hydrostatic pressure. ’ 
PRINT *,' ' 
PRINT *,’Enter Youngs Modulus, E(psi) =’ 
READ *,E 
PRINT *,’Enter Poissons Ratio, V =’ 
READ *,V 
PRINT *,’'Enter Toroid radius, R(in) =’ 
READ *,RB 
PRINT *,’Enter Circular Section radius, r(in) =' 
READ *,RL 
PRINT *,’Enter Shell thickness, h(in) =’ 
READ *,H 
PRINT *,’Enter hydrostatic pressure, P(psi) =' 
READ *,P 
PI = 4.0*ATAN(1.0) 
N= 19 
A = RB/RL 
Al=2.O*PI*E*H/(1. -V**2) 
A2=-2.Q0*PI*P*RL**2 
A3=(-PI*E*(H**3) )/(6.0*(RL*¥*2)*(1.0-V*¥*2) ) 
C GENERATE THE D COLUMN 
D(1)=A2*PI*A 
D(2)=A2*PI/2.0 
DO 30 I=3,N 
D(I)=0.0 
30 CONTINUE 
C GENERATE THE F MATRIX 
DO 40 I=1,N 
PRINT *, ‘INTEGRATING MATRIX COLUMN’ ,I 
Hell me =sAPa IMP (PAC PI, V,A,1,0.0)+A3*SIMP(FA],PI,V,A,1,0.0) 
RC Ze —=TATASIMP( FA, P1,V,A,1,1.0)+A3*SIMP(FAL,PI,V,A,1,1.0) 
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40 


F(3,1) 
F(4,1) 
ECS. 2) 
F(6,1) 
F(7,1) 
pce, 1) 
F(9,1) 


PisSiMEChAwri.V,A, 1. 2. 
A1*SIMP(FA,PI,V,A,1,3. 
SUcSIMeChASeIaV.A. I, 4. 
MLAS TMP CPA Pio VAs le 5.. 
PIXSIMECRA PIV A, 1,6: 
RIASUME (Paar, V,A,1,7- 
PEXSIMPCPAS PE, VA, 1,8: 


O)+AS3*SIMPC(PAL, FPI,V,A,1,2. 
18) 
.Q) 
0)+A3*SIMP(FA1,PI,V,A,I,5. 
0)+A3*SIMP(FA1,PI,V,A,1,6. 


0)+A3*SIMP(FA1,PI,V,A,I, 3 
0)+A3*SIMP(FA1L,PI,V,A,1I,4 


0)+A3*SIMP(FA1,PI,V,A,I,/7 
0)+A3*SIMP(FA1,PI,V,A,1,8 


0) 


0) 
0) 


.O) 
.0) 


Melo) = AL*SIMPCFA,PI,V,A,1,9. 


eri.) 
Hei 1) 
HG 3,1) 
F(14,1) 
Pilon) 
EG@ibjy) 
PEC 1.) 
FiGES, L) 
Riese t) 
CONTINUE 


A1*SIMP(FB,PI,V,A,I,1. 
A1*SIMP(FB,PI,V,A,1,2. 
A1*SIMP(FB,PI,V,A,1,3. 
A1*SIMP(FB,PI,V,A,1I,4. 
A1*SIMP(FB,PI,V,A,1,5. 
A1*SIMP(FB,PI,V,A,1,6. 
A1*SIMP(FB,PI,V,A,1,7. 
A1*SIMP(FB,PI,V,A,1I,8. 
A1*SIMP(FB,PI,V,A,1,9. 


O)+A3*SIMP(FAL,PI,V,A,1,9. 
0)+A3*SIMP(FB1,PI,V,A,I,1. 
0)+A3*SIMP(FB1,PI,V,A,I1,2. 
0)+A3*SIMP(FB1,PI,V,A,1I,3 
OUFAS*SIMP(FBL, PI(V,A,1,4: 
0)+A3*SIMP(FB1,PI,V,A,1I,5. 
0)+A3*SIMP(FB1,PI,V,A,1I,6 
OMA eSeSIMECEB1 , PI,V,A,1, 7 
OTASeSiMechol. ri, V,A,1,8. 


0)+A3*SIMP(FB1,PI,V,A,1,9. 


CONVERT THE LINEAR SYSTEM TO A TRIANGULAR SYSTEM 
CALL GAUSS(F,D,N, ERROR) 
IF (ERROR) THEN 
PRINT *,'MATRIX GENERATED IS SINGULAR, ' 
Peat e* SOLUTION iS NOT POSSIBLE.” 


ELSE 


SOLVE THE TRIANGULAR LINEAR SYSTEM 
CALEB SOLVEGE DN, C) 
PRINT *,'THE FOLLOWING ARE THE 10 COEFFICIENTS FOR THE’ 
PRINT *, ’NORMAL DISPLACEMENT OF THE TOROID. '‘ 


4-0 


, 


PRINT *,’USE IN THE SERIES; An*COS(n*theta/r) ' 


PRINT *,’ 
Pet ee 


An’ 
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An n 


0) 
0) 
0) 


10) 


0) 
0) 


=O) 
<8) 


0) 
0) 








C 


EROND GX (le, 6 (1), f=1),,10) 

PRINT *, ' ' 

PRINT *,’THE FOLLOWING ARE THE 9 COEFFICIENTS FOR THE’ 
PRINT *,’'TANGENTIAL DISPLACEMENT OF THE TOROID. ’ 

PRINT *,’ ° 

PRINT *,’USE IN THE SERIES; Bn*SIN(n*theta/r) ’ 

PRINT +, 

PRINT *,! n Bn n 


su Bn’ 


PRINT *,(1-10,€C(T), l=11,19) 
ENDIF 
END 


C THIS SUBPROGRAM FUNCTION INTEGRATES THE MATRIX FUNCTIONS 
C DEFINED BY USING SIMPSON'S RULE AS AN APPROXIMATION 


100 


REAL FUNCTION SIMP(F,PI,V,A,N,Z) 
REAL PI,F,H,SUMEVN,SUMODD,X,B 
INTEGER I 
EXTERNAL F 
B=0 .0 
H=PI/100. 
SUMEVN=0 .0 
SUMODD=F(A,V,N,Z,H) 
DO 100 I=1,49 
X=2 .*I*H 
SUMEVN=SUMEVN+F(A,V,N,Z,X) 
SUMODD=SUMODD+F(A,V,N,Z,X+H) 
CONTINUE 
SIMP=(H/3.)*(F(A,V,N,Z,B)+4. *SUMODD+2 . *SUMEVN+F(A,V,N,Z,PI)) 
RETURN 
END 
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C DEFINE THE MATRIX FUNCTIONS 
C CALCULATE ROWS 1 THROUGH 10 
C FUNCTION FA REPRESENTS MEMBRANE ENERGY 
REAL FUNCTION FA(A,V,1,2Z,X) 
REAL Q 
Teer ct. Ll) THEN 
Q=I-1 
FA=(A+COS (X) )*COS (Q*X) *COS (Z*X ) +2 . O*V*¥COS (X) *COS (Q*X) * 
+ COS (Z*X)+(COS (X) **2 ) *COS (Q*X) *COS (Z*X) / (A+COS (X) ) 
EESE 
Q=1-10 
FA=(A+COS (X) ) *Q*COS (Q*¥X) *COS (Z*X) +V¥* (COS (X) *Q*COS (Q 
+ *X) *COS (Z*X) -SIN(X) ¥SIN(Q*X) *COS (Z*X) ) -COS (X) *S1 
+ N(X) *SIN(Q*X) *COS (Z*X) /(A+COS (X) ) 
ENDIF 
END 
C FUNCTION FAl REPRESENTS BENDING ENERGY 
REAL FUNCTION FA1(A,V,1,Z,X) 
REAL Q 
Tomei el. 11) THEN 
Q=I-1 
FA1=(A+COS (X) ) *( (Q*Z) **2) *COS (Q*X) *COS (Z*X) -V*SIN(X 
= )* (Q¥ (Z**2 ) ¥SIN (Q¥X) *COS (Z*X)+ (Q¥*2) ¥Z*COS (Q¥X) * 
+ SIN(Z*X) )+(SIN(X) **2) *Q*Z¥*S IN (Q*X) *SIN(Z*X) / (A+C 
+ OS (X)) 
EUSE 
Q=1-10 
FA1=(A+COS (X) ) *Q* (Z**2) *COS (Q¥X) *COS (Z*X) - V*SIN(X)* 
+ (Q*Z*COS (Q*¥X) *SIN (Z*X) + ( Z**2) ¥SIN (Q*X) COS (Z*X) ) 
+ +(SIN(X)**2)*Z*SIN(Q*X) *SIN(Z*X) /(A+COS (X) ) 
ENDIF 
END 
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C CALCULATE ROWS 11 THROUGH 19 
C FUNCTION FB REPRESENTS MEMBRANE ENERGY 
REAL FUNCTION FB(A,V,I,Z,X) 
REAL Q 
IF (I .LT. 11) THEN 
Q=I-1 
FB=(A+COS(X) )*Z*COS (Q*X) *COS (Z*X) +V* (COS (X) *Z*COS(Q 
a *X) *COS (Z*X) - SIN(X) *COS (Q*#X) *SIN(Z*X) ) -COS(X) *SI 
- N(X)*COS (Q*X) *SIN(Z*X) / (A+COS (X) ) 
ELSE 
Q=I-10 
FB=(A+COS (X) ) *Q*Z*COS (Q#X) *COS (Z*X) - V*SIN(X) * (Z*SIN 
“ (QkX) *COS (Z*X) +Q*COS (Q#X) ¥SIN(Z*X) )+(SIN(X) **2)* 
se SIN(Q*X) *SIN(Z*X) /(A+COS (X) ) 
ENDIF 
END 
C FUNCTION FB1 REPRESENTS BENDING ENERGY 
REAL FUNCTION FB1(A,V,I,Z,X) 
REAL Q 
IF (I .LT. 11) THEN 
Q=I-1 
FB1l=(A+COS(X) )*Z*(Qk*2) *COS (Q#X) *COS (Z*X) -V*¥SIN(X)* 
+ (QkZ*SIN(Q#X) *COS (Z*X) + (Q#*2) *COS (Q*X) ¥SIN(Z*X) ) 
e +(SIN(X)**2) *Q*SIN(Q*X) *SIN(Z*X) /(A+COS (X) ) 
ELSE 
Q=I-10 
FBl=(A+COS (X) ) *Z*Q*COS (Q*X) *COS (Z*X) - V*SIN(X) *(Z*SI 
“ N(Q*X) *COS (Z*X) +Q*COS (Q#X) *SIN(Z*X) )+(SIN(X) **2) 
% *S IN (Q*X) ¥SIN(Z*X) /(A+COS (X)) 
ENDIF 
END 
C 
C GAUSSIAN ELIMINATION MODULE 
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SUBROUTINE GAUSS(A,B,M, ERROR) 
INTEGER M,I,J,K, INDEX 
REAL A(M,M) ,B(M) , PIVOT, TEMP, RATIO 
LOGICAL ERROR 
ERROR=. FALSE. 
DO 100 I=1,M 
INDEX=I 
PIVOT=ABS(A(I,I)) 
DO 200 J=I+1,M 
IF (ABS(A(J,I)) .GT. PIVOT) THEN 
PIVOT=ABS(A(J,I)) 
INDEX=J 
ENDIF 
200 CONTINUE 
IF (INDEX .GT. I) THEN 
DO 400 K=I,M 
TEMP=A(I,K) 
A(1,K)=A( INDEX, K) 
A( INDEX, K)=TEMP 
400 CONTINUE 
TEMP=B (I) 
B(I)=B( INDEX) 
B( INDEX) =TEMP 
ENDIF 
IF (PIVOT .EQ. 0.0) THEN 
ERROR = .TRUE. 
ELSE 
DO 300 J=I+1,M 
RATIO=A(J,1)/A(I,1) 
DO 500 K=I+1,M 
A(J,K)=A(J ,K)-A(I,K)*RATIO 
500 CONTINUE 
B(J)=B(J)-B(I)*RATIO 
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300 


100 


C 


CONTINUE 
ENDIF 
CONTINUE 
END 


C SUBROUTINE TO SOLVE TRIANGULAR LINEAR SYSTEM 


100 


200 


SUBROUTINE BSOLVE(A,B,M,Z) 
INTEGER M,1,J 
REAL A(M,M) ,B(M) ,Z(M) , SUM 
DO 200 I=M,1,-1 
SUM=B(T1) 
DO 100 J=I+1,M 
SUM=SUM-A(I,J)*Z(J) 
CONTINUE 
Z(1)=SUM/A(I,I) 
CONTINUE 
END 
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